    subroutine compute_station_deviations(c,nd,NewNet,MIN_WA_COMP)
!
! Programmed on Jun-17th-2013 by F.M.
!
! Compute various kinds of avearages (per station) of the deviations of the 
! station-magnitudes from event magnitudes.

    include 'DiBonaNet_structure.FI' 
    include 'seisev_EVENT_PAHSE_structure.FI'
    parameter (imaxv=6000)
    integer nd,i,j,k,MIN_WA_COMP
    real*4 v(imaxv),w(imaxv)
    real*4 rmean, standev, median
    logical ok,weighted_mean
    
    write(*,*)' subroutine: compute_station_deviations'
! I have to make all the computation station by station.
! Then, for each station, I have to pass all the catalog searcing for good
! data. The maximum number of available data per staion is about 2600. So I dimension
! the vectors at 6000.
    
    open(file='station_average_deviations.csv',status='unknown',unit=5)
    write(5,'(a)')'code,sc,scerr,n_data,mean_dev,'// &
                  'mean_dev_sd,wgt_mean_dev,'//      &
                  'wgt_mean_dev_sd,median_dev,huber_dev,huber_dev_sd'
    
    do j=1,NewNet.nsta
      k=0
      NewNet.sta(j).n_data = 0
      do i=1,nd    
        if(c(i).EVok .and. NewNet.sta(j).code .eq. c(i).sta .and. c(i).WAok)then
           k=k+1
           v(k)=c(i).MLDB_dif
           w(k)=c(i).MLDB_s_wgt
        endif
      enddo
 !
 ! memo del numero di dati a disposizione: 
      NewNet.sta(j).n_data = k
      
 ! compute MEAN 
      if(k.gt.MIN_WA_COMP)then
        call mean_standev(k,rmean,standev,v)
        NewNet.sta(j).mean_dev = rmean
        NewNet.sta(j).mean_dev_sd = standev
      else
        NewNet.sta(j).mean_dev = 0.
        NewNet.sta(j).mean_dev_sd = 0.
      endif 
      
 ! compute WEIGHTED MEAN
      if(k.gt.MIN_WA_COMP)then
        rmean=0.
        standev=0.
        ok=weighted_mean(k,rmean,standev,v,w)
        if(.not.ok)then
          rmean=0.
          standev=0.
        endif
        NewNet.sta(j).wgt_mean_dev = rmean
        NewNet.sta(j).wgt_mean_dev_sd = standev
      else
        NewNet.sta(j).wgt_mean_dev = 0.
        NewNet.sta(j).wgt_mean_dev_sd = 0.    
      endif       

 ! compute MEDIAN
      if(k.gt.MIN_WA_COMP)then
        call mediana(v,k,median)
        NewNet.sta(j).median_dev = median
      else
         NewNet.sta(j).median_dev = 0.
      endif
     
 ! compute HUBER TRIMMED MEAN    
      if(k.gt.MIN_WA_COMP)then
        ruse = 0.3
        call Huber_dev(k,v,ruse,rmean,wsd,wmax,wmin,w)
        NewNet.sta(j).huber_dev = rmean
        NewNet.sta(j).huber_dev_sd = wsd
      else
        NewNet.sta(j).huber_dev = 0.
        NewNet.sta(j).huber_dev_sd = 0.
      endif
      
      if(k.gt.MIN_WA_COMP)then
        WRITE(5,'(a5,a1,2(f7.3,a1),i5,a1,7(f7.3,a1))')   &
          NewNet.sta(j).code ,',',    &
          NewNet.sta(j).sc ,',',    &
          NewNet.sta(j).scerr ,',',    &
          NewNet.sta(j).n_data ,',',    &
          NewNet.sta(j).mean_dev ,',',    &
          NewNet.sta(j).mean_dev_sd ,',',    &
          NewNet.sta(j).wgt_mean_dev ,',',    &
          NewNet.sta(j).wgt_mean_dev_sd ,',',    &
          NewNet.sta(j).median_dev ,',',    &
          NewNet.sta(j).huber_dev ,',',    &    
          NewNet.sta(j).huber_dev_sd
       endif
  
    enddo
    
    close(5)
    
    return
    end